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1 Introduction 



Driven dissipative systems arise in modelling many different phenomena in Nature, 
e.g. avalanches, forest fires, earthquakes and Darwinian evolution. Such systems 
frequently exhibit long range temporal and spatial correlations and scaling laws 
analogous to statistical mechanical systems at a critical point. The concept of self- 
organized criticality was introduced [0J to describe systems where critical behaviour 
arises without the fine tuning of any parameter. 

The study of self-organized criticality is essentially a study of a family of dy- 
namical systems which evolve in time in small or large steps which have a power 
law distribution and we shall call "avalanches" but they may have a different inter- 
pretation. The canonical example of a self-organized critical system is the abelian 
sandpile model Jl|, || where an integer variable Zj is assigned to each site i in a lattice 
and the dynamics is given by choosing a random site i in the lattice, increasing Zj 
by 1, Zi — > Zi + 1, and if zi exceeds a threshold value z c , then an avalanche start and 
the value of Zi is decreased by redistributing some of Z, among the neighbours. This 
in turn may push some of the Zj over threshold, where j is a neighbour of i, and so 
on. These rules are then supplemented by appropriate boundary conditions. See 
for a recent detailed discussion. 

It is natural to introduce a mean field theory for self-organized criticality M. 
This has been done in many different ways for different models yielding identical 
critical exponents || ||, 0, ||. One of the approaches is to study sandpiles on a Bethe 
lattice ||. In this case different parts of an avalanche front are uncorrelated and 
one obtains the scaling law 

N w (k) ~ k-^ 2 (1) 

where N w (k) is the frequency of avalanches involving k sites. We shall refer to k as 
the width of the valanche. 

In this paper we introduce a simple model for the propagation of a one-dimensional 
avalanche front. This model was suggested by studying earthquakes in the Burridge- 
Knopoff model || |K| [11, 12), The basic idea is to assume that neighbouring 



parts of an avalanche are correlated in the following way. If we label elements of the 
avalanche by integers % and an element i moves a distance hi then the displacement 
of its neighbour, labelled by i + 1, is distributed with a probability distribution 
Pi(h i+ i) = 4>{h i+ i — hi) which is centered on hi but otherwise independent of i. This 
distribution is modified by an appropriate boundary condition at h i+ i = amount- 
ing to a certain killing probability for the avalanche. In this paper we shall in fact 
work with the simplest possible probability distribution cf) which corresponds to the 
avalanche front performing a Bernoulli random walk on the positive integers and 



2 



terminating once it returns to 0. In this case one can perform explicit calculations 
and we find an exponent 3/2 in the scaling law relating the frequency of avalanches 
to their width. 

If we define the size of an avalanche as 

A = 5>. (2) 

i 

then we find that the frequency of avalanches of size A is given by 

N(A) ~ A' 4 / 3 . (3) 

Sometimes, e.g. in earthquake models, it is natural to place an upper bound on 
the maximal allowed displacement hi and in that case the power law turns into an 
exponential decay for large A. The relation to earthquakes is discussed in more 
detail in [O . 



2 The random walk model 

In this section we define the random walk model. We introduce a class of random 
walks (or paths) and each path will correspond to an avalanche of a particular 
form. Questions about the size, shape and frequency of avalanches can therefore 
be translated into questions about the number of paths satisfying the appropriate 
conditions. In this model we are not able to address questions about the correlations 
between different avalanches nor discuss the self-organized critical spatial structure 
that arises in realistic models. 

We shall refer to the discrete elements participating in an avalanche as blocks. 
Consider a semi-infinite chain of blocks Mi, lying along the x-axis, labelled by the 
non-negative integers which can be regarded as the x coordinates of the blocks. We 
assume that the blocks move in integer steps in the direction of the y axis with 
the x coordinate unchanged. In applications to many one-dimensional systems, e.g. 
earthquakes and sandpiles, the x coordinate would be the x coordinate of the blocks 
at time while the y coordinate, in our present terminology, would be the addition 
to the x coordinate in an avalanche. 

We imagine that the block labelled by i = 1 is the first one to move in an 
avalanche and the next block to move is M2 etc. For convenience we assume that 
the block at the origin Mq is fixed at all times, but this assumption is not at all 
essential. We could easily have another avalanche moving to the left and this second 
avalanche would be uncorrelated to the one moving to the right. We assume that Mi 
moves a distance 1 along the y-axis. This is the initial condition for an avalanche. 
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We now assume that M 2 moves a distance 2 along the y-axis with probability ^ and 
a distance with probability |. In general, given that the block M n has moved 
a distance h n , block M n+ i moves a distance h n ± 1 with equal probabilities. By 
definition, the avalanche stops once the next block to move in fact does not move. 
The minimal avalanche is the one where M 2 does not move and this is the most 
likely event with probability |. It is sometimes natural to place an upper bound h 
on how far the blocks can move so we modify the rules by requiring M n+ \ to move a 
distance ft, — 1, if M n has moved a distance h. Of course it would be natural to allow 
neighbouring blocks to move the same distance. This would not affect any of the 
conclusions but make some of our equations more complicated. In this formulation 
an avalanche can be regarded as a directed random walk on the non-negative integers 
which is reflected from h and stops when it returns to 0. 
We can write 

n 

h n = Yl a ii ( 4 ) 
i=l 

where G\ = 1 and <7j for % > 2 is a random variable which takes the values ±1 with 
equal probabilities. The duration T of the walk is defined by 

h T = 0, h n >0, 0<n< T. (5) 

The size (moment in seismology) A of the avalanche is given by the area under the 
graph of the function h n , i.e. 

A=Y, K. (6) 

n=i 

We are interested in calculating the frequency of avalanches with size A and the 
relation between A and the width T. 

3 Short time behaviour 

We begin by considering the model in the absence of an upper cutoff h, i.e. we put 
h = oo and return to the case of finite h later. A directed random walk (or path) is 
one where the x coordinate increases by 1 in each step. Let W denote the set of all 
directed walks in the positive quadrant of the xy-plane which start from (0, 0) and 
return to the x-axis. Such a walk is located at (1, 1) after one step. Let N(A,T) 
denote the number of paths in W which return to the x-axis after T steps and whose 
graph, together with the x axis, encloses an area A given by Eq. (||), see Fig. 1. Let 

N(T)=Y,N(A,T). (7) 

A 
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The probability that the walk returns to the origin after T steps is given by 



P(T) = 2~ T+1 N(T). (8) 

We divide by 2 T_1 since the first step in the walk is given. Similarly, the conditional 
probability P(A\T) that a walk covers an area A, given that it lasts a time T, can 
be written as 

It follows that the probability that an avalanche has size A is given by 

P{A) = 2- T+1 J2N(A,T). (10) 

T 

It is well known from the solution of the classical gambler's ruin problem, see e.g. 
(T§, that 

P(T) = — x [ T Tj l ) (11) 
if T is even and otherwise. For large T 

P{T)~T~i. (12) 

Since we are discussing random walks it is natural to expect the average height (h n ) 
of a walk which lasts a time T to scale like VT for large T. The average area A 
should therefore grow as T 3 / 2 for such walks. A priori one would expect entropic 
repulsion to play a role here so 3/2 should be a lower bound on the "average area 
exponent" but, as demonstrated below, there is in fact no shift away from the naive 
value of this exponent. 

If we consider directed random walks that are allowed to cross the x axis and 
define the area under the walks to be positive if the walk is in the upper half plane 
and negative when it is in the lower half plane, we can use ([|) and (^) to express the 
area as a sum of independent but non-identical random variables. The generalized 
central limit theorem |16[ applies to this sum and we find that asymptotically the 
area is normally distributed around zero with a variance T 3 . Assuming that P(A\T) 
is normally distributed around its average value with a variance T 3 we expect that 
for large A we have 

P(A) = J2 p (A\T)P(T) 

~ Jo T exp [—r* — ) dT 

~ A~ 4/3 . (13) 
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Below we shall verify that the asymptotic behaviour of P(A) is indeed given by Eq. 
(f[3|) even though one can prove that P(A\T) is in fact not normally distributed by 
computing its first few moments. 

Let W denote the class of directed walks in W which avoid the line y — 1 until 
they return to y — 0, i.e. if w G W and if returns at time T then w;(x) > 1 for 
1 < x < T — 1, where denotes the ^-coordinate of the path for x = t. Let us 
denote by N(A, T) the number of paths in W which return to at time T and cover 
an area A. Then 

N(A,T) = N(A-T + l,T-2), (14) 

see Fig. 2. Now consider any directed walk w G W which lasts a time T > 2 and 
covers an area A. Let 7\ denote the smallest integer > 1 such that w(T\) = 1. The 
largest possible value of T\ is of course T\ = T — 1. If we cut the path u> in two 
pieces at the point (T 1; 1) then we can associate uniquely to w two paths, w G W and 
u>i G W, of duration T% + 1 and T — Ti + 1, respectively, see Fig. 3. In the extreme 
case T\ — T — 1 the second walk is the trivial one of length 2. If we denote the area 
under the first walk by A\ then the area under the second one equals A — A\ + 1 
and we find that 

T-l A 

iV(A T) = 5 A1 5 T2 + E E T\ + 1)N(A - A x + 1, T - T x + 1) 

Ti=lAi=l 
T-l A 

= <W5t 2 + E E ^(Ai-T^Ti-l^A-Ai + ^T-Ti + l), (15) 

Ti=lAi=I 

by Eq. (|14] ) . The first term on the right side of Eq. (|l^) corresponds to the two step 
path. We define the generating function f(z,u) for the numbers N(A,T) by 

oo oo 

/(^l^^iV^T)//, (16) 

T=2A=\ 

which is convergent for \z\ < 1 and \u\ < \. Eq. ( p!5| ) can now be rewritten as 

f(z, u) = zu 2 + f(z, u)f(z, uz). (17) 

For z = 1 we can easily solve this equation and find 

1 1 



/(i )U ) = r -vr^, (is) 

which for u = ~ takes the value | in accordance with Eq. ([10|) . 

For general values of z and u Eq. ( pT7| ) is not explicitly soluble. We note however 
that the equation can be rearranged to read 

'<*■»> = T=7&ST (19) 
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which, upon iteration, yields a continued fraction expansion! for f(z,u) 

/(*,«) = 3"^- (20) 



1 - 



z 3 u 2 



Looking at Eq. 
we expect 



we expect the average size of an avalanche to diverge, i.e. 



•i dz 



oo. 



2 = 1 



Indeed, we shall prove that 



Qn 

dz n 



f(z,u) 



(21) 



(22) 



2 = 1 



as u | 9, for any n > 1. Let us define 



P u (A)=^ M T - 1 iV(A,T). 

T 



Then 



2 = 1 



For any path of duration T, covering an area A, it is easy to see that 

-T — 2<A< -T 2 . 
2 ~ ~ 4 



(23) 



(24) 



(25) 



The lower bound is obtained by considering the path which zigzags between 1 and 
2 and covers the smallest possible area while the upper bound corresponds to the 
triangular path that climbs to height T/2 in time T/2 and then descends to zero in 
time T/2. It follows that for u smaller than but close to | we have the bounds 



e - Cl ^- 2u)A P(A) < P U (A) < e - c *- 2 ")^P(A) 



(26) 



4 It is perhaps of interest to note that the continued fraction (g0|) appeared in a letter from 
Ramanujan to Hardy written in 1913 |17 where it was used to express some remarkable identities, 



one of which can be written, in our notation, 



'5 + VE VE+i 



,277/5 



However, the Ramanujan identities seem unrelated to the statistical properties of the random walk 
model we are interested in. 
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where c\ and c 2 are positive constants. It is natural to regard u as a temperature- 
like parameter and critical point so we expect scaling as this point is 
approached. Assuming that 



P U (A) ~ - 4u 2 A /3 ) (27) 

for large A, where F is a function decaying more rapidly than any power, we find 
that 

J2A n P u (A) ~ (1 - 4« 2 )T. (28) 

A 

Since this is valid for any n it follows from Eq. fl2"2"|) that 



/3=i and 7 = |. (29) 

This is of course consistent with the inequalities fl26D which imply that | < /3 < ~. 

In order to verify Eq. (p2]) let us denote the derivatives of the generating function 
/ with respect to the first and second argument by d\f and d 2 f, respectively. We 
begin by considering the case n — 1. Differentiating Eq. fll7D with respect to z we 
obtain 

d\f(z,u) = u 2 + f(z,uz)d 1 f(z,u) + f(z,u)(d 1 f(z,uz) + ud 2 f(z,uz)). (30) 

Putting z = 1, rearranging and using Eq. we obtain 

u 2 + ud 2 f(l,u) 



dif{l,u) 



l-2/(l,u) 
u 2 + 2u 2 (l-4u 2 H 



1 (31) 



1 -4m 2 



Assume now that Eq. ( p2| ) holds for n < iV — 1 where N > 2. Differentiating Eq. 
(|17j) N times with respect to z yields 

^/(l,«) = E(fc )#/(!> «) £ ( y»- i - k dldp- i - k f(l,u). (32) 

Rearranging we find that 

d?f(l,u) = l - Y ( ) &f(l,u) 

l-2/(l,u) t'o \ k I 



X 



N-k , N _ k 



E \ U^-^^- fe /(l^). (33) 



J 
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By the inductive hypothesis the most singular terms on the right hand side of Eq. 
(P^) correspond to j '• + k = N if k > and j = N — 1 when k = 0. The desired 
result follows. 

Working slightly harder we can determine the coefficient of the leading divergence 
of the Nth moment of P U (A) and this allows us to place a further restriction on the 
function F introduced above. In view of Eq. fl22|) one can write 



d?f(l, u) = C N (1 - An 1 )*'- + 0((1 - 4U 2 ) 1 --). (34) 



It is straightforward to check that Eq. (^) determines the following recursion rela- 
tion for the coefficients Cn 

N-l 

k 



C n = E ( ^ J CkC N - k + (3N - 4)iYC , JV _ 1 . (35) 



fe=i 

It follows that up to power corrections 

C N ~ (2N)\ (36) 

for large N. 

Suppose now that the function F in (2?) is an exponential function of a power, 

i.e. 

F{x) = e~ axq (37) 

for some constants a,q > 0. Using the ansatz ([Z7]) to calculate the iVth moment of 
the area distribution we find that q is fixed to equal 3/2 by the asymptotic formula 
(|36l). We therefore expect that 

P U {A) ~ A-^e-^ 1 - 4 " 2 ) 37 ^ 172 ( 38) 

in agreement with the bound (Effi). 

This completes our disussion of the frequency- size distribution of avalanches 
where we do not need to take the upper cutoff h into account, i.e. the exponent 4/3 
governs the size distribution of small avalanches in the presence of a cutoff. We now 
turn to the study of the tail of the size distribution and will see that for large A 
the probability P(A) falls off exponentially when the upper bound h is taken into 
account. We also estimate the area for which we have a transition from a power law 
to an exponential decay. 



4 Long time behaviour 

We now consider a directed random walk with a reflecting barrier at height y = h. 
Let Pi(t) be the probability of the walk being at height i after t steps. Then the 
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initial condition is 

Pi(l) = 5 a . 

We let p(t) denote the column vector whose ith entry is Pi(t). Then 



p(t) 



-M 



t-i 



p{\) 



(39) 



(40) 



where M is the matrix 



M 



( 


1 













o\ 


2 





1 
















1 





1 
















1 





1 

















1 














1 








V o 











1 


o / 



(41) 



The probability that a walk lasts exactly T steps is evidently 

P{T)= l -p 1 {T-l). 

Let (•, •) denote the standard inner product on and e^, i = 
orthonormal basis. Then we can write 



P(T) = (e , ( -M 



T-l 



ex). 



Let D denote the matrix whose elements are defined by 



D, 



z l 5. 



(42) 

h, the standard 
(43) 

(44) 



i,j — 0,...,h. If P(A,T) denotes the probability that a walk lasts a time T and 
covers an area A then P(A, T) is given by the coefficient of z A in the matrix element 

(eo^MDf-'e,). (45) 

The generating function for the probabilities P(A, T) can therefore be expressed as 

u 2 MD 



Q(z,u) = (e , 



(46) 



1 - uMD 

since the Neumann series for the inverse of (1 — uMD) is easily seen to converge for 
\u\ < \ and \z\ < 1. 

Eq. (|^) allows us in principle to calculate P(A,T). However, the interesting 
feature of P(A, T) is that it falls exponentially with T and P(A, T) = unless 
A < hT. It follows that P(A) = J2t P{A,T) falls exponentially with A provided A 



10 



is large enough. In order to establish this exponential decay it suffices to show that 
Q(l,u) is finite for some u > |. We can write 

g(1 ' M) = ^ (e °'A^M ei) (47) 

where A = (2m)" 1 . Evaluating the matrix element in Eq. (f47|) is an elementary 
calculation and we find 

Q{1 , u)= l^m (48) 

4 cos ((h — 1)0) 

where 

e 



16 - A + iv/T^A 1 (49) 



and we are assuming A < 1. The first singularity of Q(l,u) as u moves beyond | 
is encountered for the smallest 9 G [0, 27r) for which the denominator in Eq. (f48l) 
vanishes, i.e. 

" =2^' (50) 
It follows that the radius of convergence of Q(l, u) is 



r = lf^W^)' (51) 

and for large A we find 

P(A) < Ce~ cA/h2 (52) 

where c and C are positive constants and c can be taken to be independent of h. 

The exponential decay of P{A) takes over from the power law found in the 
previous section for A « h 3 since a random walk must have at least h 2 steps in 
order to feel the effect of the reflecting barrier at y — h. In order to prove this note 
that we can write 

P{T) = f ^T dU > (53) 
where the contour encloses the unit disc in the complex plane. Calculating the 
residues we find that 



P(T) = 2,r> Y S in^ ^' +2 "> cos- ^ + 2 "> (54) 

n=0 fl 1 fi 1 

and consideration of this formula at large h shows that P(T) crosses over from 
exponential decay to a T~ 3//2 decay for T ^ h 2 . 
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5 Discussion 



By universality, we do not expect the principal results we have obtained to change if 
we replace the simple Bernoulli random walk, considered in this paper, by a random 
walk with any rapidly decaying transition function. The exponent value 4/3 ought 
to be universal as will be exponential decay for large A in the presence of an upper 
cutoff h. 

We do, however, expect to be able to change the power law decay by considering 
strongly correlated random walks. In earthquake models, for example, the size 
distribution exponent for small and intermediate size events varies from around 2/3 
to values greater than 1. It is clear that by looking at the statistics of avalanches in 
a realistic system one can concoct a random walk model with an identical avalanche 
distribution. 

The more interesting problem of understanding correlations between different 
avalanches cannot be studied in the random walk framework unless one introduces 
different interacting random walks. The principal virtue of the model we have dis- 
cussed is that it gives us a qualitative and quantitative insight into the genesis of 
power law distribution for avalanches without introducing any complicated dynam- 
ics. 
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Fig. 1 A directed random walk of duration T = 12. 
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Fig. 2. This figure illustrates the one to one correspondence between paths in W 
of duration T and paths in W of duration T — 2. 
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Fig. 3. This figure illustrates how one can uniquely decompose any directed path 
into a pair of paths in W' and W. 
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